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We provide a systematic test of empirical theories of cova- 
lent bonding in solids using an exact procedure to invert ab 
initio cohesive energy curves. By considering multiple struc- 
tures of the same material, it is possible for the first time 
to test competing angular functions, expose inconsistencies 
in the basic assumption of a cluster expansion, and extract 
general features of covalent bonding. We test our methods 
on silicon, and provide the direct evidence that the Tersoff- 
type bond order formalism correctly describes coordination 
dependence. For bond-bending forces, we obtain skewed an- 
gular functions that favor small angles, unlike existing mod- 
els. As a proof-of-principle demonstration, we derive a Si 
interatomic potential which exhibits comparable accuracy to 
existing models. 

PACS numbers: 61.50.Lt,34.20.Cf,33.15.Dj 

Large-scale atomistic simulations are becoming in- 
creasingly important in the study of complex physical 
phenomena such as fracture, plastic deformation, two- 
and three-dimensional melting, epitaxial growth, shock 
wave propagation, friction, sintering, etc. Ideally, one 
would like to represent the atomic interactions in these 
simulations with a quantum mechanical approach, treat- 
ing explicitly the electronic degrees of freedom. This is 
a computationally demanding proposition, tractable at 
present only for relatively small system sizes of order 10^ 
atoms. An alternative description is in terms of effec- 
tive interatomic potentials which allow fast evaluation of 
energies and forces, making possible simulations involv- 
ing more than 10* atoms. The drawback in going from 
an explicit quantum treatment of electrons to an effec- 
tive interatomic potential is a significant loss in accuracy 
that may undermine simulation results. In most cases, 
the microscopic mechanisms of greatest interest (for in- 
stance, bond formation or rupture) are precisely those 
which require high degree of transferability, that is, abil- 
ity of the potential to describe accurately a wide range 
of local atomic environments. 

Over a decade of experience has shown that such 
transferability is difficult to attain, especially in covalent 
solids, for inherently quantum effects such as bond bend- 
ing and breaking, hybridization, charge transfer, and 
metalization. In the prototypical case of silicon, about 
thirty model potentials exist in the literature (l) , includ- 
ing popular and innovative ones by Stillinger and We- 
ber (SW) §, Tersoff @, and Chelikowsky et. al. [|. 
Although the shortcomings of existing model potentials 
have been carefully documented, it has proven very dif- 
ficult to improve them or to understand, even qualita- 



tively, the causes of their failures [|l],^ . Some theoretical 
arguments have been advanced to motivate the form of 
an effective interaction [^,^ and to derive potentials as 
approximations of quantum models p|-pT| , but little spe- 
cific theoretical guidance exists to aid in the development 
of new potentials. The most successful approach to date 
is to guess a functional form using physical intuition and 
then adjust parameters to fit a database of ab initio struc- 
tural energies 0]. The reliance on intuition and fitting 
leads to the two questions that motivate our work: (1) 
Is there an ab initio justification for the functional form 
of an interatomic potential? and (2) Given a particular 
form, is there a systematic way to obtain new potentials 
directly from ab initio energy data? 

In this Letter, we present an exact procedure for in- 
verting ab initio energy data to obtain parameter-free 
many-body potentials JT^ . The inversion approach was 
pioneered by Carlsson, Gelatt, and Ehrcnrcich (CGE) for 
the case of a pair potential JlSf , and since then the same 
formula has been applied with limited success by only 
a couple of authors @Jl^. We revisit the inversion ap- 
proach with the following innovations: (i) a recursive for- 
mulation that incorporates many-body interactions and 
strains other than uniform volume expansion; (ii) the re- 
quirement that ab initio energies be exactly reproduced 
for relevant densities only (near the equilibrium solid and 
liquid densities); and (iii) the use of an overdetermined 
set of structures for the same material, which guaran- 
tees a wide range of relative atomic arrangements. These 
ideas form a general framework for analyzing functional 
forms and deriving potentials, as illustrated by applica- 
tion to Si. In this manner, we provide satisfactory an- 
swers to both questions posed in the previous paragraph. 
We analyze the two defining features of covalent bonding 
in two steps, 1. Pair bonding and 2. Angular forces. 

1. Pair bonding: We begin with the simplest case of 
a pair potential, in which the cohesive energy E[(j)\ of an 
arbitrary structure is given by. 



i^j p=i 



(1) 



with atomic separations grouped into shells Sp of radius 
Spr containing Up atoms each. Dilation of the lattice is 
achieved by varying the parameter r with the structural 
quantities {sp} and {n^} fixed. Shells are numbered so 
that si < S2 < S3 < . . ., and distances scaled so that 
si = 1. A simple rearrangement of the terms in Eq. (Q) 
yields the desired inversion formula for 
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0(r) = - Eir) - Y: 

1 \ p=2 



np(l){spr) 



(2) 



Although the unknown potential appears on both sides 
of this equation, recursive substitution generates the ex- 
plicit formula. 
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(3) 



which was originally derived by CGE invoking the linear- 
ity of the functional E[(j)] [^. Our recursive formulation 
generalizes to nonlinear functionals and suggests a simple 
computational procedure: If the tail of 0(r) is assumed 
known for r > a, then the potential is uniquely deter- 
mined by solving the recursion in order of decreasing r 
starting at r = a (because Sp > I for p > 2). An impor- 
tant case is that of finite range, i.e. 0(r) = for r > a, 
as is typically assumed for interatomic potentials. 




FIG. 1. The inverted pair potential for silicon (i) before and 
(ii) after a cutoff is imposed, compared with ipswir) (dashed 
line). Numbers inside the figure indicate shell radii in the 
diamond lattice. The inset shows the diamond LDA data and 
the interpolant (i) before and (ii) after imposing a cutoff. 

To illustrate the inversion procedure, we apply it to an 
ab initio database consisting of cohesive energy curves 
for Si from density functional calculations in the local 
density approximation (LDA) fl^ ]. In order to keep the 
procedure simple while still capturing the important lo- 
cal bonding characteristics, the database includes the fol- 
lowing crystals: (i) the low-energy and low-coordination 
structures, diamond (Si-I), ^-tin (Si-II) H, BC-8 (Si- 
III) [0, and BCT-5 |l|; (ii) SC and FCC crystals for 
metallic behavior; and (iii) the graphitic structure for 
non-tetrahedral hybridization 1 19 1 . These structures have 
coordinations 4,6,4,5,6,12, and 3, respectively. We con- 
sider only atomic volumes smaller than (3.54A)'^ to avoid 
the difficulty of LDA to represent accurately the energies 
of isolated atoms pO[ . Smooth interpolation of the LDA 



data and extrapolation to infinite volume with an ex- 
ponential tail are used. The LDA data points for the 
diamond lattice with the interpolant are shown in the 
inset of Fig. 1. 

The inverted pair potential for the diamond curve, 
shown in Fig. 1, is clearly unphysical: Its long range and 
strong repulsion at the first neighbor distance contradict 
our intuitive understanding of covalent bonding. Similar 
' results have been obtained in previous work applying the 
CGE formula to metals 0,|l^ and semiconductors [ pT[ . 
Our recursive approach reveals that these problems are 
inherent to the inversion process, which, in spite of being 
exact, stretches the assumption of a volume-independent 
potential to an unphysical extreme. Because inversion 
amounts to solving in order of decreasing distance from 
infinite separation, the tail of the potential comes from 
unscreened interactions between atoms in a low density 
gaseous phase . The same tail is then used to describe 
long range interactions in a bulk crystal, which are pre- 
sumably screened by the presence of closer atoms. While 
the nature of screening and its description by an effec- 
tive potential are subjects of active research, it is obvious 
that distant atoms in a bulk crystal cannot interact in the 
same way as atoms with the same separation in a gas. 




FIG. 2. Inverted pair potentials (with cutoff) for seven sil- 
icon bulk phases. The inset shows the implied bond order 
p extracted from these curves (points) compared to y/4jz 
(line). p{l) reflects the 812 bond length and energy Q. 

To rectify the inversion procedure, we forgo the re- 
quirement that the potential exactly reproduce the entire 
cohesive energy curve. Instead, we focus on condensed 
volumes typical of solid and liquid environments, whose 
exact energies can be preserved with any choice of tail for 
the potential. For Si, we find that exponential decay to a 
cutoff near the second-neighbor distance in the diamond 
lattice, 3.84A, generates potentials in good agreement 
with bonding theory. For example, as shown in Fig. 1, we 
force the energy curve to be zero for r > asw = 3.77118A 
without disturbing energies within 10% of the equilib- 
rium bond length, where covalent bonds are well-defined, 
to produce an inverted pair potential with a deep mini- 
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mum at the first neighbor distance. 

Applying the same procedure to the other curves in 
our database, we obtain the potentials of Fig. 2. The 
large discrepancy between them is direct evidence for a 
well-known fact: the energetics of silicon cannot be de- 
scribed by a pair potential alone 0. Our results suggest, 
however, that an environment-dependent pair potential 
can describe the ideal bulk phases reasonably well. There 
is a clear coordination dependence to the curves: bond 
lengths (positions of the minima) increase, and bond 
strengths (depths of the minima) decrease with increas- 
ing coordination. This behavior can be described by the 
bond-order formalism, which is justified on grounds of 
theoretical arguments ^-[gj] as well as experience with em- 
pirical potentials |]3|,^,^3[. In its simplest form, a bond 
order potential is given by. 



^W = EE E 9iRMRjk)h{e,jk), (5) 



(r, Z)^Mr)+p{Z)4>A{r), 



(4) 



where 4>r and (pA are monotonic repulsive and attractive 
terms, respectively, and p{Z) gives the bond strength as 
a function of the coordination Z. The leading order ap- 
proximation of the bond order is p{Z) oc Z~^^^. Since 
this comes from describing the local density of states 
by the bandwidth only, we expect the approximation to 
work well for the metallic phases with Z > A (BCT-5, (3- 
tin, SC, and FCC). For the covalent phases with Z < 4 
(diamond, BC-8, and graphite), band shape effects be- 
come important and we expect significant departure from 
the Z~^/'^ behavior. 

If the repulsive interaction (pR were known, the bond- 
order term could be extracted directly from the ab initio 
data, using p{Z) = VA(ro)/V"/"(ro), where Va = (j)-(j)R, 
ro is the minimum of the inverted potential cf), and we 
set p = 1 for the diamond lattice {Z = 4). The repulsive 
term, intended to represent an effective force between 
electrons due to Pauli exclusion, is the weakest link in 
bond-order models, since its form must be assumed and 
then fit to empirical data without theoretical guidance. 
Although the general trend is insensitive to the choice 
of we find that using cfyR = , where cf)^ is 

the repulsive part of the SW potential, produces a p{Z) 
which lies remarkably close to its expected behavior (see 
inset of Fig. 2) . 

2. Angular forces: While the energetics of bulk phases 
can be fairly well described by a bond order pair func- 
tional, it is well-known that many-body interactions with 
explicit angular dependence are required for silicon, for 
example, to stabilize the diamond lattice against shear 
strain ||,^. As the simplest case of a many-body po- 
tential, we consider one with volume-independent pair 
terms and separable three-body terms like the potentials 
of SW and Kaxiras and Pandey Q . The many-body en- 
ergy, F{r) = E{r) — V2(r), formed by subtracting the 
pair terms V2{r) from the total energy E{r), is expressed 
as a sum over pairs of bonds. 



where cos 9ijk = Rij ■ Rjk- Following theory ||^,|7| and 
practice we assume F > 0, which implies that the 
pair potential come from the diamond lattice inversion 
described above. A particular form for the angular term 
h(9) must be assumed in order to invert F[g, h] for the 
radial function g[F,h]. The procedure is the same as in 
the pair potential case: solve Eq. (||) for g(r) to obtain a 
recursion. Grouping bonds into shells as above and tak- 
ing the positive root of the resulting quadratic equation 
yields the desired expression, 



-/3(r) + V/3(r)2+4a„(f(r)-7(0) 
9{r) = ^ , (6) 



2aii 



where 

pq — J2r,jl£S; 



a. 



and 7(r) = Y.'^^2J2'^=p^pq9{spr)g{sqr). In the app 
sums, only k > j contributes to avoid double counting. 
An explicit formula like the CGE pair potential can be 
obtained by recursive substitution with Eq. (^) and in- 
volves a tree of nested square roots. It is simpler to follow 
the same computational procedure as before, solving the 
recursion in order of decreasing distance starting at the 
cutoff. In principle, the same general procedure can be 
applied to determine radial functions of other forms or for 
higher order terms in a cluster expansion of the effective 
potential. For example, for a nonseparable three-body 
term involving three bond lengths, like the potential of 
Pearson et al. |Q, the recursion comes from solving a 
cubic equation, and for a four-body term, a quartic equa- 
tion. 




FIG. 3. An inverted angular function for silicon from 
the BC-8, BCT-5, and /3-tin energy curves compared with 
hsw(d)- The inset shows the collapse of the inverted g{r) with 
the average curve (dashed line) and gswif) (dotted line). 

As in the pair potential case, it is useful to invert 
more than one cohesive energy curve of the same ma- 
terial for g[F,h]{r). If the assumed angular dependence 
h{9) (and two-body terms) were truly transferable, then 
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the same radial function g{r) would result from every in- 
version. Conversely, the greater the variance between the 
inverted g{r), the less transferable is the assumed h{9). 
This principle gives us a quantitative means of assess- 
ing the quality of angular functions directly from the ab 
initio data. For example, although the SW angular func- 
tion h{6) = (cos(0) -I- |)^, produces mediocre collapse of 
the inverted g{r), it is fortuitously far better than taking 
h{9) — {cos{9) + i)*, due to the latter's flatness near the 
tetrahedral minimum. 

Using the same principle, we can extract an optimal 
angular dependence from the ab initio data by assum- 
ing a series expansion, h{0) — J2^^QCi{cos{9) + |)^^'- 
For the curve shown in Fig. 3 (defined by cg = l,ci = 
— 1.86,C2 = 1-42), the collapse of radial functions from 
the low energy phases is rather good, as seen in the in- 
set of Fig. 3. A novel feature of the inverted angular 
function is its skew about the minimum to favor smaller 
angles. This is consistent with the conclusion that exist- 
ing potentials tend to overpenalize angles smaller than 
7r/2 [|], which presumably leads to poor descriptions of 
surfaces, clusters, and certain defects. The skewed an- 
gular function also raises the energy of overcoordinated 
metallic structures relative to covalent ones by penalizing 
large angles. While it is typical to characterize metallic 
structures by the presence of small angles Q], we note 
that metallic structures tend to have angles near tt also. 
Covalent bonds are actually characterized by angles in 
the intermediate range, 7r/2 to 27r/3. 

Although we are not attempting here to provide an im- 
proved potential for Si, we have performed some tests of 
the potential obtained by the inversion just described p^ ] 
(in this proof-of-principle demonstration, we omit coor- 
dination dependence for practical reasons). We find that 
it performs as well as the popular SW and Tersoff poten- 
tials without fitting to any defect structures, for energies 
of other silicon bulk phases, defects such as interstitials 
and vacancies, generalized stacking faults, the concerted 
exchange diffusion mechanism and (100) and (111) 
surface reconstructions. Considering the database em- 
ployed in the inversion, we conclude that many important 
features of chemical bonding are contained in cohesive 
energy curves for ideal bulk phases. 

In conclusion, we have presented a general procedure 
for inverting cohesive energy curves to obtain many-body 
effective interatomic potentials. By inverting ab initio 
cohesive energy curves for silicon, we have demonstrated 
how general features of bonding are revealed. Elsewhere 
we will describe extensions of these ideas, for example, to 
the inversion of energy curves for shear strains to obtain 
the angular function h[F,g](9) directly. The inversion 
procedure provides a systematic method for deriving in- 
teratomic potentials and a unique tool for understanding 
their general limitations through the direct use of ab ini- 
tio data. It is hoped that this tool will lead to potentials 
with improved transferability, a goal that has proven elu- 



sive when pursued by intuitive arguments and fitting of 
databases. 
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